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EVALUATION OF A NEW METHOD OF INTEGRATING THE ORBITAL 

EQUATIONS OF MOTION FOR USE IN SPACE NAVIGATION 

By John D. McLean 
Ames Research Center 


SUMMARY 


The applicability of a new method of integrating the orbital equations of 
motion, developed by Dr. J. M. A. Danby of Yale University Observatory, to the 
problem of space navigation is investigated. The investigation is carried out 
by means of a digital computer program written for the method. The transearth 
and translunar phases of a circumlunar trajectory are taken as sample problems, 
but the method can also be applied to any orbital mission. The gravitational 
effects of the sun, moon, earth, and earth's oblateness are considered. 

A detailed description of the computer program, including listings, is 
presented. Data show the effects of various error sources and the trade off 
between computing speed and accuracy. The method is compared with Cowell's 
method on the basis of computer storage, integration time, and accuracy. 


INTRODUCTION 


In recent years considerable study has been devoted to the problems of 
navigation and guidance for manned space missions. Except for near-earth sat- 
ellites, the astronaut will probably be provided with the capability, at least 
as a secondary system, for carrying out these tasks without the aid of ground- 
based equipment. If such an on-board navigation and guidance system is not to 
be restricted to limited emergency procedures, a digital computer will be 
required aboard the spacecraft. 

One of the more demanding functions of such a computer is to solve the 
vehicle's equations of motion. The most commonly proposed method is the numer- 
ical integration of the differential equations of motion, although there are 
other possibilities such as the interpolation of stored data. Two well-known 
methods of setting up the equations to be integrated are Cowell's method in 
which the total accelerations acting on the spacecraft are integrated, and 
Encke's method in which differential equations for perturbations from an oscu- 
lating conic are solved. Both of these methods require a complex integration 
routine, and any alternative which will substantially reduce the time and stor- 
age requirements is desirable. For on-board use the reduction in storage 
requirements is particularly desirable because the decrease in the number of 
components enhances reliability and reduces weight and power requirements. In 
addition, the problem of accurately starting the numerical solution of a system 
of differential equations is quite complicated. Since the estimated trajectory 
for a navigation system must be restarted every time new observational data are 
obtained, it would be desirable to avoid or minimize this difficulty. 



The purpose of this report is to present an evaluation, from the stand- 
point of application to on-board computation, of a new approach to the problem 
of integrating the equations of motion. This integration procedure was devel- 
oped by Dr. J. M. A. Danby of Yale University Observatory and is described in 
detail in references 1 and 2. As in Encke’s method, the perturbations of the 
true orbit from a reference conic are computed, but the method reduces the inte- 
gration of the perturbing accelerations from the solution of a set of differen- 
tial equations to simple quadratures. Although a well-known mathematical 
technique is used to obtain the perturbations, Danby was, to the author f s 
knowledge, the first to recognize the advantage of applying this technique to 
orbital equations. This method of solving the equations of motion will be 
referred to in the remainder of the report as "Danby f s method." It will be 
shown how Danby T s method can be used to reduce the problems discussed in the 
preceding paragraph. 

The original application of Danby f s method, as described in references 1 
and 2, fits the trajectory with a small number of conics, called mean orbits. 
These mean orbits closely approximate the true orbit, including perturbations, 
thereby allowing it to be studied using closed form equations. The advantages 
of this procedure for many theoretical studies are obvious, but the accurate 
mean orbits are unnecessary for the specialized application of space navigation. 
Since the calculation of the mean orbits requires each portion of the trajec- 
tory to be integrated several times, it is desirable to avoid such calculations. 
Likewise, the method, as described in reference 2, uses eccentric anomaly as the 
independent variable in order to avoid inverting Kepler’s equation. Since 
space navigation requires frequent processing of data at accurately known times, 
it would be desirable to use time as the independent variable. In this study 
several modifications have been made to Danby’s method, as presented in the ref- 
erences, in order to eliminate these two difficulties without sacrificing com- 
puter storage requirements and speed. These modifications are described in 
detail and the modified system is compared with Cowell’s method 1 for translunar 
and transearth trajectories. 


NOTATION 


a semimajor axis of conic 

Ay lower bound on Aq 

Aq measure of validity of Simpson’s rule 
Ap upper bound on position perturbation 

Aq lower bound on position perturbation 

Ay upper bound on Aq 

throughout this report "Cowell’s method" will refer to Cowell’s method of 
setting up the orbital equations of motion for integration rather than Cowell’s 
numerical integration method. 
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parameters of closed form conic equations 

increment in time 

increment in eccentric anomaly 

a 3x3 unit matrix 

position vector 

|r| 

vector of small deviation in position 

|r| 

total number of integration steps 
time 

6x1 matrix of perturbing accelerations 
gravitational potential 

potential of central body to which conic is referred 
velocity vector (equivalent to R) 

|V| 

vector of small deviation in velocity 

M 

Cartesian components of vehicle T s position 

6x1 matrix (state vector of position and velocity deviations 
from conic) 

transformed forcing function 
gradient operator 
increment in time 

magnitude of terminal position and velocity deviations of 
Danby solution from Cowell solution 

increment in eccentric (or hyperbolic) anomaly 
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0 


state transition matrix 


cpi,cp 2 ,cp 3 ,cp 4 submatrices of 0 


Superscripts 

T transpose of a matrix 

. derivative with respect to time 

_ a vector or 3x1 matrix 

Subscripts 

i,n integers 

o initial value, usually referring to time of last rectification 

THE INTEGRATION METHOD 


The basic principles of Danby T s method as presented in reference 1 are 
summarized here for the convenience of the reader. A more general coverage of 
the mathematical principles involved can be found in reference 3 or in texts 
dealing with vector and matrix differential equations or multidimensional con- 
trol theory. 

The basis of the method is the solution of the orbital equations of motion 
in terms of a reference conic and a set of associated perturbations . The com- 
putation of the perturbations is accomplished by simple quadratures through 
the use of the state transition matrix, in contrast to Encke’s method which 
requires the solution of a system of differential equations. The transition 
matrix, $(t, t Q ) , is the matrix of first partial derivatives of the components 
(all the trajectory calculations discussed in this report are carried out in 
Cartesian coordinates) of position and velocity at time t with respect to 
the same quantities at time t Q . It can be shown that if x is a 6x1 matrix, 
or state vector, of position and velocity deviations from the reference trajec- 
tory then 

x(t) = ®(t, t Q )x(t 0 ) (1) 

The use of this matrix in the integration is explained as follows: 

The orbital equations of motion may be written as 


R + F = VU(R + r ) + u(R + r, t) (2) 

where U is the gravitational potential due to the central body of the refer- 
ence conic, R is the position vector on the conic, r is the vector of 


k 



position deviations of the true orbit from the conic and u is the vector of 
perturbing accelerations. If it is assumed that r is sufficiently small, 
the true orbit may be approximated in terms of linear (first order) perturba- 
tions from the reference conic. In this case 

F - F x (R)r + u(R, t) (3) 

where Fi is the 3x3 matrix of first partial derivatives of VU with respect 
to the components of R. 

Since R is a function only of position on the reference conic, hence of 
time, equation ( 3 ) represents a system of three linear second-order differen- 
tial equations which can be rewritten 

x = 

where 


u 


and 

F = 

It can be shown that the solution of the homogeneous part of equation (4) is 
given by equation (l) and that the complete solution is 


x - <&(t, t 0 )x 0 + $(t, t Q ) J t 0 )u( T )dr (5) 

t o 

Since the reference trajectory is a conic, the elements of 0 can be 
computed in closed form. 2 Furthermore, u is a function only of the position 
vector on the conic and time; that is, it does not depend explicitly on the 
perturbations. Thus the integrand in equation (5) is a known function of time 
and a simple quadrature formula, such as Simpson 1 s rule, can be used to solve 
for each component of x separately. 

2 Since 0 also satisfies the homogeneous part of equation (4); that is, 

i = F$ 

where $(t 0 ) = I, 0 may be found by numerical integration. In this case the 
reference orbit need not be a conic and the effects of perturbing accelera- 
tions on the transition matrices may be included. This is the method used to 
obtain 0 in Cowell program B discussed later in the report. 


as a first-order system to give 

Fx + u (4) 
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In addition, the problems Involved in starting or restarting a numerical 
Integration of a differential equation are eliminated. The elimination of a 
complex restart procedure is particularly Important in a navigation system 
because the integration is stopped and restarted frequently to allow processing 
of observed data and revision of the estimated trajectory. 

As the integration proceeds, the perturbations will grow until r is so 
large that equation (3) is no longer valid. When this occurs it is necessary 
to obtain a new reference conic (a process called rectification) for which 
equation (3) is valid. It is customary in most integration schemes to rectify 
only after r becomes too large, in which case the Integration must be 
restarted at a point where r is still below the desired bound. It is pos- 
sible, however, to rectify any time before r exceeds its bound. 


APPLICATION TO ON-BOARD COMPUTATION 


The purpose of the present study is to evaluate the utility of Danby T s 
method for use in on-board space navigation systems, and a digital computer 
program was written for that evaluation. The details of the program, which 
was designed to be incorporated Into a complete (simulated) guidance and navi- 
gation system at a later date, are presented in the appendixes. Appendix A 
contains the formulas used for the reference conic and a derivation of the 
equations (different from the equations of refs. 1 and 2) for the transition 
matrices. The equations for the components of the forcing function and the 
quadrature formula are given in appendixes B and C, respectively, and the pro- 
gram listings are given in appendix D. 

As was stated earlier, certain modifications were made in the mechaniza- 
tion in order to use Danby’s method to best advantage In the space navigation 
problem. These modifications are discussed in detail in the following para- 
graphs . 


The Reference Orbit and Rectification 

The first modification Is in the choice of the reference conic and this 
in turn requires a different rectification procedure. In the application 
described in references 1 and 2, the integration is started at time t G using 
the osculating conic as the reference orbit and rectification is carried out 
at some time, t x , when r exceeds a predetermined value. At that time a new 
reference conic Is computed which has a different initial velocity and passes 
through the position on the perturbed orbit at t x . The integration and rec- 
tification are repeated several times, always starting at t Q and going to 
progressively later values of ti, until a reference conic, or mean orbit, is 
obtained which approximates a large segment of the true orbit very closely (in 
position) . Because of this repeated integration, it is unnecessary to con- 
strain r to very small values, and a large number of integration steps are 
made between rectifications. 
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Since the determination of an accurate mean orbit is not necessary for 
space navigation, it was decided to eliminate the repeated integration and use 
the osculating conic as the only reference orbit. In this case a more strin- 
gent restriction must be put on r in order to obtain the desired accuracy, 
and it was found that over most of the trajectory rectification must occur 
after each integration interval. If the value of r is used to determine 
when to rectify, the maximum allowable value must be exceeded before it is 
known whether rectification is necessary. Thus the new conic must originate 
at the beginning of the last integration step, and the integration over this 
last step must be repeated. If this rectification procedure is used with the 
osculating conic as the reference trajectory, most of the trajectory will be 
integrated twice. For this reason the program was arranged to require recti- 
fication after each integration step. There is negligible penalty in computing 
time or storage for such frequent rectification, but it was found, . as will be 
discussed later, that excessive round-off errors are encountered near the cen- 
ters of attraction. This difficulty can be avoided while still maintaining 
efficient use of computing time by the use of double precision for the computa- 
tion of the reference conic. An alternate method would be to take several 
integration steps between rectifications near the centers of attraction while 
rectifying after every integration step elsewhere. 


Choice of the Independent Variable 

The second modification in the use of Danby’s method was in the choice of 
the independent variable. Danby (ref. 2) recommended that the eccentric anom- 
aly of the Keplerian orbit be used in order to avoid the iteration necessary 
for the inversion of Kepler’s equation. While this choice has obvious advan- 
tages for a trajectory study, it is not particularly suitable for an on-board 
navigation system. Such a system requires that quantities from which the tra- 
jectory can be determined be observed at intervals. By some method a trajec- 
tory is estimated and the values that the observations would have if the space- 
craft were on the estimated trajectory are computed. Because of observational 
errors the actual trajectory can never be known, and that estimated trajectory 
which minimizes the differences or residuals between the actual and computed 
observations is found. Since the times of the observations are fixed, they 
will occur at different eccentric anomalies for each estimated trajectory. 
Hence, the eccentric anomalies at the times of the observations are unknowns 
to be determined, while the times associated with them can be known quite 
accurately . 

This difficulty presumably could be overcome by some sort of interpola- 
tion scheme, but the cost in computer storage would be large. In addition, if 
a Kalman filter type (ref. 4 ) of trajectory determination is used, the transi- 
tion matrices between observations would be needed, and these must relate the 
states at the two different times. For these reasons, time was chosen as the 
independent variable. However, in order to minimize the number of solutions 
of Kepler’s equation, the integration was set up as outlined below. 
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The Integration 


Although the program in its final form was set up to rectify after each 
integration step, this discussion is presented so as to be equally applicable 
to the case where several integration steps are made between rectifications. 
(The program presented in appendix D can be modified fairly easily to take 
several steps between rectifications.) As an aid to clarity in discussion, t Q 
is defined as the time of the last rectification while 0 is the change in 
eccentric anomaly between t Q and a subsequent time. 

Assume that equation ( 5 ) has been integrated from t Q to t n and it is 
desired to continue the integration to t n + 2 . This integration proceeds as 
follows : 

(1) At t n , the elements of u(t n ) and $(t n ,t 0 ) have been computed in 
the previous step unless t n = t 0 . In the latter case 0 is the unit matrix 
and u(t n ) must be computed. 

( 2 ) The forcing function u is transformed to time t Q by 

y n " 0 *k 0 )u(t n ) 

where the transformed forcing function y is the integrand in equation (5)* 

(3) The angle 0 is increased from 0 n to (0 n + hjg) and the associated 
position and velocity are computed along with the time, t n+1 , the transition 
matrix o(t n+1 , t Q ) , and u(t n+1 ). Then 

y n+1 = t 0 )u(t n +i) 

( 4 ) The angle 0 is increased to (0 + 2 hg) and the conic is extended to 
t n+2 . If ^n+2 is smaller than (or equal to) the next time at which an obser- 
vation, velocity correction, or termination of the flight is to occur 

y n +2 = t n+s) u ( t n+2) 

is computed. This amounts to making two equal increments in 6 and accepting 
the associated change in time. If t n + 2 is greater than the next desired time 
of stopping, the value of (0 n+2 - 0ti) which will make the two increments equal 
is found by iteration. This change in eccentric anomaly is halved to give a 
new hg and the process is repeated from (2). (in the remainder of the report 
the term "step size" will refer to hjj; and "integration step" to the interval 
< t < t n+2 . ) 

(5) The y n , which constitute a set of values of the forcing function at 
times t n , t n+1 , and t n+2 transformed to the rectification time, t Q , are 
integrated by quadratures. 

These are sufficient data for the quadrature of each component of the 
forcing function by Simpson’s rule to give 
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®~ X (r, t n )u(T)dT = a oy n + a!y n+1 + a 2 y n +2 (6) 


but the quadrature is complicated by the fact that the time increments are 
unequal. Simpson 1 s rule is derived on the assumption that the curves to be 
integrated can be fit with a quadratic polynomial in the independent variable 
which, in this case, is time. If the time intervals are equal, the coeffi- 
cients of the polynomial are constants, but for unequal intervals they must be 
calculated each time. Expressions for these coefficients are derived in 
appendix C. 

Adjustment of the Step Size 

The incremental eccentric anomaly h^ over which the quadratic approxi- 
mation is valid changes markedly over the length of a lunar trajectory. This 
change is most radical near the centers of attraction and can be attributed 
mainly to the rapid change of the elements of the transition matrices in this 
region (see ref. 5)- Thus, for greatest efficiency it is desirable to adjust 
the step size as the integration proceeds. If the eccentric anomaly were the 
independent variable, the fourth differences of the y n could be computed and 
would give a good indication of the accuracy of the quadrature (see ref. 6). 
The step size could then be adjusted to be the maximum size which would 
restrict the magnitude of the fourth difference, and hence, of the error, to 
the desired limit. 

It is possible to derive an expression which is equivalent to the fourth 
difference for unequal time increments, but the resulting computation is 
rather cumbersome. Since simplicity is one of the objects of this investiga- 
tion, it was decided to try a simpler procedure which was found empirically to 
be quite satisfactory. This procedure is explained with the aid of sketch (a) 



Sketch (a) 
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The area under the quadratic curve in the sketch represents the ith 
component of the integral given by equation (10) . Since the perturbing accel- 
erations are known to be smoothly varying functions of time, it was reasoned 
that the amount of curvature in the quadratic curve passing through the y n ± 
could be used as a measure of the size of higher order fluctuations. To assume 
y n to be a constant would be the simplest approximation and would result in 
an^ integral equal to the area C^. A linear approximation would add the area 
to this integral while the quadratic curve would add the area given by the 
algebraic sum (Aj[ + Bq) . A measure of the curvature of the quadratic curve is 
therefore the ratio of Ai to (Aj_ + B±) . The six components of the integral 
were considered collectively as follows: Let rq 2 be the sum of the squares 

of the Aj_ representing position deviations while Vq 2 is the sum of the 
squares of the Aj_ representing velocity deviations. Similarly, let r ^ and 
vijt 2 represent the sums of squares of (Aj_ + B±) • Finally, define Aq as 



The step size is automatically adjusted so that the Aq remains in the range 

a l < a q < A u 

where Ay and Ay are input constants. 

The step size must be small enough for the linearity assumption (that the 
true orbit can be expressed accurately in terms of small perturbations from 
the reference conic) to still be valid. This means that the error, Sr, in the 
position perturbation vector, r, must be small compared to the position vector 
of the reference conic. If one uses the error analysis of reference 1 on a 
simple one -dimensional example, it is seen that over most of the trajectory 
Sr will be roughly proportional to r 2 /R 2 , where R 2 is the distance from 
the vehicle to the perturbing body. However, when the distance, R x , to the 
central body is small compared to R 2 (particularly when the earth is the cen- 
tral body), Sr may become proportional to r^/Rj. . Since it is desired to 
bound Sr/Rj. , which is proportional either to r 2 /R 2 2 or r 2 /R x R 2 , it was 
decided for simplicity to require 


A s < I < a r 

where R is the smaller of R x and R 2 , and A3 and Ap> are input constants. 3 

The program in appendix D is arranged so that if either Aq or r/R 
exceeds its upper bound after a given integration step, the step size is halved 
and the integration step is repeated. If both quantities drop below their 
lower bounds, the step size is doubled for the next interval. Part of the 
digital computer study discussed in the next section was devoted to the problem 
of choosing values for the parameters Ag, Ap>, Ay, and Ay which control the 
value of hy . 

3 If it is not desired to rectify every integration step, then a good 
criterion is to rectify when r/R exceeds Ar. 
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TEE DIGITAL COMPUTER STUDY 


This section of the report presents the results of a digital computer 
study used to evaluate the application of Dauby* s method outlined in the pre- 
vious section. The first part of the discussion deals with the choice of 
desirable values for the parameters Ag, Ag, and Ajj and some of the data used 
in establishing these values are used to illustrate the effects of round-off 
error. It was found that a r must be chosen on the basis of a trade-off 
between speed and accuracy, and the second part of the discussion deals with 
this trade-off and with the time required to invert Kepler *s equation. 

Finally, Danby*s method is compared with Cowell's method in terms of the time 
and computer storage required to perform various tasks. 


Choice of System Parameters 


The first step in the numerical study was to determine desirable values 
for the upper and lower bounds, Ag, Ar, Ar, and A\j discussed earlier. These 
bounds establish a band of allowable error for each integration step, but the 
location of the actual error within this band is random. For example, the 
position deviation, r, from the reference conic and the corresponding error 
in r increase nonlinearly with time and the rate of increase also varies 
along the trajectory. Sketch (b) illustrates how this can affect the accuracy 
of the integration. Suppose the nth integration step covers the time from 


r/R 



Sketch (b) 

t Q to t3 so that the perturbation (shown in the sketch by the solid line) 
from the reference conic, and hence the error due to nonlinearity, nearly 
attains its maximum allowable value. Now consider a change in Ar. This 
might affect the previous part of the trajectory so that the nth integration 
step covers the time interval from t x to t 4 (dashed line) in which case r/R 
exceeds Ar, and hg is reduced. With hg reduced the nth integration 
step terminates at t 2 while the (n+l)th step (dotted curve) covers the time 
from t 2 to ts. With the reduced value of hg the perturbation and cor- 
responding error due to nonlinearity remain much smaller than when the nth 
integration step begins at t Q . On the other hand, while reducing hg 
reduces the error due to nonlinearity, it increases the possibility of round- 
off error. When the step size and corresponding perturbations are made quite 
small, then part of the perturbation is smaller than the least significant 
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figure in the corresponding component of the conic position or velocity and 
is lost as round-off error during rectification. (in fact, the perturbation 
could be made so small that the conic position and velocity -would not be 
altered by rectification.) 

In the interest of efficient use of computing time it would be desirable 
to make the width of the error band zero and set the upper bounds to give the 
maximum error consistent with the desired accuracy. If one attempts to 
approach this situation by the use of upper and lower bounds which are nearly 
equal, the program will back up an excessive number of times for step size 
reduction. On the other hand, if the upper and lower bounds are too far apart 
h E will be kept unnecessarily small. Thus it was decided to choose the lower 
bound which would give the most efficient operation for each corresponding 
upper bound and to use the latter parameters to establish the accuracy of 
integration. The following paragraphs describe the method used in establishing 
Ag, A p, and Ay and show the effect of round-off error on accuracy. 

Choice of lower bounds .- To establish a reasonable value for Ag, a number 
of transearth trajectories were computed with the step size adjusted only by 
r/R. (Reasonably accurate results could be obtained by proper choice of the 
initial step size.) Figure 1 shows Sp, which is defined as the total number 
of integration steps including repeats for the reduction of hp, plotted as a 
function of Ag/Ap for various values of Ap . The value of Sp is propor- 
tional to the integration time and thus is a measure of the efficiency of the 
program for a given value of Ap; Sp is not extremely sensitive to Ag/Ap 
for the larger values of Ap, but as that parameter is reduced, a definite 
minimum is obtained for a ratio of about l/5* As a result of these data, Ag 
was constrained to be l/5 of the value of Ap for the remainder of the study. 

With Ag/Ap fixed, it is now possible to determine a similar relation- 
ship between Ap and Ag. In figure 2 Sp is plotted for the transearth tra- 
jectory as a function of Ap/Ag for various values of Ap and Ag. No optimum 
value of Ap/Ag can be seen from these data, but it is clear that the smaller 
ratios require an excessive number of steps. Likewise, most of the curves 
flatten out or rise slightly for the largest ratio. These results imply that 
Ap should be about l/lO the value of Ag. 

Choice of Ay. - Because of the randomness of the errors in integration 
discussed previously, one integration is not sufficient to establish the accu- 
racy associated with a given set of values of Ap and Ag and it is desirable 
to obtain some estimate of the maximum likely error. In the trajectories used 
for figure 2, it was found that for given values of Ag and Ap, the greatest 
differences in terminal position and velocity usually occurred between tra- 
jectories having the maximum and minimum values of Ap. For this reason the 
only values of Ap/Au used in determining Ag were 0.01 and 0.25* The val- 
ues of Ap and Ag were the same as those for figure 2, except for one addi- 
tional value (0.01) for Ag. This set of trajectories was integrated with the 
osculating conic computed in double precision because when single precision is 
used the influence of Ag is obscured by round-off errors. The X component 
AX (see appendix B for the definition of the coordinate system) of the termi- 
nal position deviation from the Cowell solution for these sets of trajectories 
is plotted in figure 3 as a function of Ag for different values of Ap. The 
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X component comprises most of the deviation for this particular trajectory 
and was used to preserve sign information. The results for the single- 
precision reference conic are also presented for use later in the discussion 
of round-off error. 

Since the use of double precision essentially eliminates round-off errors , 
the errors in these solutions must be attributed to nonlinearity and inaccu- 
racies arising from the use of Simpson 1 s rule. For Ap < 10 ~ 4 the deviations 
of the solutions for Danby T s method (double precision) from the Cowell solu- 
tion are always less than 0.5 km. The question of a "correct" solution when 
the deviations are so small is somewhat nebulous and neither method can be 
considered as more accurate than the other. A considerably more detailed 
study of both Cowell T s and Danby's methods would be required to resolve the 
accuracy question further, and this study was not made since these results are 
satisfactory for navigation and guidance. 

The lack of an absolutely correct answer does not prevent the use of 
these data in establishing the effects of various error sources and determin- 
ing a reasonable value for %• For the two smallest values of A R it is 
evident that reducing Ay from 0.2 to 0.1 improves the accuracy of the solu- 
tion, but further reduction has negligible effects. The same trend is also 
seen in the curves for A R = 10 4 , but it is obscured by increased spread 
between the solutions for maximum and minimum Ap that results from nonlin- 
earity. On the basis of these results it was decided that the values of Ap 
should be set at 0.1 for best accuracy and efficiency. 

Effects of round-off errors.- We now compare the results of the single 
precision computations (of the reference conic) with those of using double 
precision in order to assess the effects of round-off error. The double and 
single precision solutions corresponding to the same combinations of Ap, Ay, 
and Ap used identical sequences of hg over the entire trajectory. For 
this reason the pertur bat ions computed at each step for the two solutions are 
nearly the same, and the terminal deviation of the single precision solution 
from that for double precision must be attributed mainly to round-off. 

The differences between these curves are within the accuracy limits of 
the Cowell integration for most of the trajectories, but a definite inaccuracy 
due to round-off is evident for all single precision solutions with A^=2xl0“ 5 
(fig. 3(a))- On the other hand, when Ap is increased to 2xl0“ 4 (fig- 3(d)), 
the single precision results are more accurate than those for double precision. 
This indicates that for Ap = 2xl0“ 4 the errors in the perturbations resulting 
from nonlinearity have become comparable to the least significant figure of 
the single precision conic position or velocity. In this case the part of the 
perturbation lost in round-off is entirely erroneous so that the round-off 
"errors" sometimes improve the accuracy. 

Consequently, it appears that by the proper choice of Ap and Ap, it is 
possible to obtain good accuracy without resorting to double precision. How- 
ever, in the vicinity of the centers of attraction, the constraint on Aq 
reduces r/R well below the value of Ap because of the rapid change in the 
transition matrices mentioned earlier. The round-off errors caused by this 
reduction in step size are evident in the curves for figure 3 when Ap = 0.01. 
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This round-off error has no significant effect in the transearth case (if 
Au>0 .01) . For the translunar case, however, the round-off errors for 
A U -0.1 are large enough to cause sizable errors. 

This phase of the digital computer study resulted in the choice of con- 
stant values of 0.2 for A^/Ap, 0.1 for Ap/Ap, and 0.1 for Ay. The next data 
to be presented -will show how, with these parameters fixed, it is possible to 
establish a trade-off between speed and accuracy by changing the value of Ap- 


Accuracy and Integration Times 

Using the values of Ag, Ap, and Ap determined in the previous section, 
it is now possible to consider the trade-off between integration time and accu- 
racy. It should be pointed out that the integration times to be presented 
here were obtained from a 709^- library clock subroutine which has a minimum 
measurable time increment of 0.6 second. 

Accuracy and integration time for trans earth c ase.- The integration times 
for the transearth trajectory are plotted 6n figure b as a function of Ap for 
both single and double precision solutions. The position deviations of these 
solutions from the Cowell results are also presented to give an indication of 
the trade-off between accuracy and integration time. 4 The velocity deviations 
in m/sec are nearly the same as the position deviations in km and have been 
omitted. Note 'that very good results are obtained for Ap < 10 ” 4 and the 
integration time is about 5-5 sec for single precision and 6.5 sec for double 
precision. Some time can be saved at the expense of errors in the order of 1 
or 2 km (and m/sec) by raising Ap as high as 5x10 “ 4 . It is doubtful whether 
one would wish to use higher values of Ap because of the rapid rise in error 
for a very small decrease in integration time. 

Accuracy and integration tim e for translunar case.- Figure 5 presents the 
integration time and position and velocity deviations for the translunar tra- 
jectories using the same values of Ap, Ap, and Ap as for the data of fig- 
ure b. The error in the single precision solution due to round-off is greater 
than 5 km for all values of a r- The results for the double precision case 
are essentially the same as for the transearth trajectory, but the deviation 
from the Cowell solution for small values of Ap is slightly larger for the 
translunar case. This result probably arises from the initial velocity being 
used differently in the two programs. Since the initial velocity is much 
larger for translunar injection, the last binary digit represents a velocity 
sufficient to produce position deviations of the order of 0. 5 km at the moon. 
The corresponding deviation for the transearth case is about l/b this amount. 


4The interaction of various parameters and initial conditions sometimes 
causes the cancellation of errors. Such cancellation produced unrealistically 
small errors for the two largest values of Ap in the single precision case 
and for Ap - 5x10 ” 4 In the double precision case. For this reason the posi- 
tion deviations given in figure A for the single precision cases are the max- 
imum encountered in the data for figures 2 and 3 with Ap set to 0.1, while 
the one point for double precision is omitted. 
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It was found that the larger errors in the single precision case arise 
because near the earth the upper bound on Aq produces very small step sizes , 
hence, very small perturbations which are mostly lost during rectification. 
This difficulty could be remedied by taking several integration steps between 
rectifications, but it was pointed out earlier that this procedure if done 
over the entire trajectory requires excessive integration times. Single pre- 
cision could be used efficiently if the program were made to take several 
integration steps between rectifications near the centers of attraction and 
to rectify every integration step elsewhere. The logic and computations for 
this procedure would require more additional storage than the use of double 
precision, but there might be a small saving in integration time. 

Time requirements for iterative solution of Kepler’s equation . - It was 
pointed out earlier that one of the big advantages of Dahby T s method over 
others for space navigation is the ease of restarting the integration. If the 
integration is to be stopped precisely at a given time in order to process 
observational data, an additional computation, namely, the inversion of 
Kepler’s equation by iteration, is required. The computation time required 
for this iteration will be different for each individual solution, but some 
representative values were obtained as follows: The initial conditions for 

the translunar trajectory were used as a starting point, and the initial 
increment of eccentric anomaly was set at 0.1 radian. The program was 
required to find the change in eccentric anomaly corresponding to a time 
increment of one hour to an accuracy of one part in 10 7 (i.e., to O.OOO 36 sec 
in this case) . The time for 100 solutions was measured, and the process was 
repeated for transearth injection. The same data were also obtained to an 
accuracy of one part in 10 s and all the results were presented in table I. 

The data for the higher accuracy were computed using double precision because 
when dt/d 0 is small the accuracy of one part in 10 s cannot be achieved with 
single precision. The accuracy actually needed remains to be established, but 
one part in 10 7 should be ample for most applications. Even an additional 
order of magnitude in accuracy does not require an excessive amount of compu- 
tation time. 


Comparison With Cowell’s Method 


Two Cowell programs were used for comparison with Danby’s method. Cowell 
program A computes only a single trajectory (no transition matrices) while 
Cowell program B computes the identical reference trajectory plus the transi- 
tion matrices. All programming for Danby’s method, except for subroutine for 
multiplying matrices, was done in the Fortran IV computer language but could 
probably be done more efficiently in assembly language (MAP) , commonly called 
machine language. On the other hand, the Cowell integration subroutine, which 
comprises nearly half of Cowell program A and over l/4 of Cowell program B, is 
written in machine language- 

It should be pointed out that the navigation and guidance task requires 
the integration to be stopped frequently for the processing of observed data 
and for the computation of several velocity corrections. These operations 
require the use of the transition matrices so that it is logical to compare 
Cowell program B, including an appropriate number of restarts, with Danby’s 
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method. However, it is of interest to know the requirements of various opera- 
tions so that data for comparison of the performance of the following tasks 
have been obtained. 

(1) The integration of a single trajectory without intermediate restarts. 

(2) The computation of a trajectory and its associated transition 

matrices without restarts . 

( 3 ) The computation of a trajectory and its associated transition 

matrices with restarts at appropriate times to simulate a naviga- 
tion problem. 

The problem of computing velocity corrections is quite complicated and 
has been left for a future study except for the following cursory examination: 
Danby's method inherently includes computing of the two-body transition 
matrices. These matrices are known to be accurate enough for use in trajectory 
determination, but introduce large errors in velocity corrections computed at 
great distances from the terminal point. There are several possible ways of 
correcting for this inaccuracy, including the use of iteration. Iteration 
would sacrifice some of the speed advantage of Danby *s method over Cowell’s, 
but would eliminate the need for a guidance reference trajectory. (Such a 
reference trajectory, discussed in reference 7; is correct from injection to 
the terminal point and should not be confused with the reference conics of 
Danby 1 s or Encke’s method.) Thus, some of the speed of the Danby integration 
would be traded for a reduction in storage requirements. Other possibilities 
include the computation of the transition matrices by the method of refer- 
ence 5; compensation for the effects of perturbations on the transition 
matrices (ref. 8), or the use of precomputed two-body aim points. 

Time requirements.- Cowell program A requires about l8 sec, exclusive of 
computer output time, to integrate the transearth trajectory and 17 sec in the 
translunar case or about 2.5 times the requirement for equivalent accuracy 
using Danby T s method. Cowell program B requires about A 3 sec for the same 
integration, that is, about 6.5 times more than is required by Danby T s method. 
The time requirements resulting from a large number of restarts will be dis- 
cussed later. 

Storage requirements .- The program using Danby’s method and single preci- 
sion throughout requires a total of 2775 words of computer storage of which 
2398 words are program and 377 are data. These figures do not include the 
Fortran monitor system or the subroutines for obtaining the ephemerides of the 
sun and moon. The ephemerides computation, which is the same for both methods 
of integration, requires about 1700 words of storage- Since it is known that 
this requirement can be reduced by about an order of magnitude (by use of a 
more specialized method), it has been eliminated from comparison. When the 
osculating conic is computed in double precision, the storage requirement 
increases by 72 data words for a total of 2849* 

For comparison, Cowell program A requires 2568 words of storage of which 
237 are for data, while Cowell program B requires 4172 words, of which 8l4 
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are for data. Thus, when used as part of a navigation system, Danby’s method 
requires about 1/3 less storage than the Cowell program. 

Effects of restarting .- Cowell program B and the single precision version 
of Danby 1 s method were modified slightly to allow the integration to be 
stopped and restarted. The transearth trajectory was then integrated with 
each program, the integration being stopped and restarted at 4-5 times which 
was considered sufficient for making observations for trajectory estimation. 
The resulting change in terminal position was 0.02 km for Cowell’s method and 
0-34- km for Danby’s method (using A p> =■ 10“ 4 , Ay = 0.1, and Ap =■ 0.01) . The 
integration time for the Cowell program increased to about 2.4- min, 5 while 
that for Danby 's method increased to about 8 sec. 

The data in table I indicate that 4-5 additional solutions of Kepler’s 
equation contribute no more than about 1 sec of the additional integration 
time for Danby’s method. This figure is confirmed when it is noted that the 
total number of integration steps, including repeats for reducing step size 
or stopping at the desired time, is increased from 83 to 126 . The total num- 
ber of integration steps could probably be reduced by the use of additional 
program logic. However, the integration time is still small, and it is doubt- 
ful whether this complication is worthwhile unless a considerably greater num- 
ber of observations is to be made. 

Furt her improvements for use in space navigation systems . - For the sample 
navigation problem just considered, it has been found that Cowell’s method 
requires about 17-5 times as much computing time as Danby’s method. The 
Cowell program can be modified to reduce this factor to between 8 and 10 at 
the sacrifice of accuracy. Likewise, the integration time for Danby’s method 
could probably be reduced further by programming improvements including the 
use of machine language. However, in the particular problem considered, it is 
known that the integration time required by Danby’s method is no greater than 
the time required to process the data from 45 observations. Further investi- 
gation is needed to establish whether the complete navigation system can best 
be improved by a more refined integration method or by more efficient data 
processing. 


CONCLUDING REMARKS 


The data presented here show that for integrating a single reference 
trajectory and computing the associated transition matrices Danby’s method 
requires only about two-thirds the storage required by Cowell’s (program B) 
method and about one-sixth the computing time. If the transition matrices are 
omitted, the storage requirement for Cowell’s method is slightly less than for 

5 0ther Ames programs using the Cowell integration for simulation of 
guidance and navigation systems (e-g., see ref. 7) use fixed step mode between 
those observations which are scheduled at reasonably short intervals. These 
programs require less integration time than the above figure, but the stops 
for observations result in terminal errors of about 6 km for this number of 
observations . 
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Danby's method, but the integration time is still about three times that 
required for Danby's method. The accuracies of the two methods are (as far as 
could be ascertained) comparable and are satisfactory for space navigation. 

The introduction of a number of stops and restarts to allow for estimation 
of the trajectory from observed data increases the time required by both 
methods. However, the increase for Cowell's method is so great that in the 
sample problem considered, the Cowell integration required 17-5 times as much 
computing time as Danby's method. The restarts have little effect on the 
accuracy of either solution, but it is possible to reduce the integration time 
for Cowell's method by a factor of about 2 at the expense of fairly large ter- 
minal errors . 

A suitable Encke program was not available for comparison. However, 
Encke's method requires the solution of the same conic equations, while the 
computation of the perturbations is more complex than for Danby's method. 

Also, since Encke's method requires the solution of a set of differential 
equations for the perturbations, a complex procedure is required to restart the 
the integration. On this basis, Danby's method can be assumed to be at least 
competitive with Encke's, particularly if the transition matrices are needed. 

Finally, while the equations for the conic solutions and the transition 
matrices are rather complex, they require only a knowledge of the calculus and 
elementary matrix theory for their formulation. Furthermore, no detailed 
knowledge of numerical methods is required to set up the quadrature of perturb- 
ing accelerations. For these reasons, Danby's method offers a simplicity and 
versatility which seems particularly well suited to engineering applications. 


Ames Research Center 

National Aeronautics and Space Administration 
Moffett Field, Calif., June 21, 1966 
125 - 17 - 05 - 01-21 
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APPENDIX A 


EQUATIONS FOR THE CONIC TRAJECTORY AND TRANSITION MATRICES 


Laplace's f and g formulation, which is summarized below, was used to 
calculate the reference conic. (See ref. 9 for a detailed explanation.) Fol- 
lowing the discussion of the reference conic, the derivation of the equations 
for the transition matrix is presented. This derivation is intended to cor- 
respond as closely as possible with the program listing (SUBROUTINE PARTLS) 
given in appendix D and for this reason may seem rather awkward from a mathe- 
matical point of view. Reference 10 presents different formulations of the 
equations for the conic trajectory and the transition matrix which are valid 
for all conics, while the ones given here break down in the parabolic case. 

The formulas in reference 10 may also require less computer time and storage 
and their use should be considered in future applications. 


The equations for the reference conic_are based on the assumption that 
the position and velocity vectors, R Q and V Q , respectively, at time t Q are 
given. Then at time t 


R =• fR Q + gV Q 

v - m , o + gv Q 


(Al) 


The following set of formulas (eqs. (A2) through (A9) ) for the scalars f, g, 
f, and g and appropriate version of Kepler's equation are given by Pines 1 in 
an unpublished work. 



1 Pines, Samuel: Analytic Mechanics Associates, Uniondale, New York. 
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where \i is the product of the universal gravitation constant and mass of the 
central body. The semimajor axis, a, is given by 


i = 2_ _ 
a R 0 “ (j. 

This equation is valid for both elliptic and hyperbolic orbits and gives a 
negative value for a in the hyperbolic case. 


The definitions above are valid for both ellipses and hyperbolas , pro- 
vided fj_ are defined as follows: 


(l) For elliptic orbits 


(2) For hyperbolic orbits 


fi = 0 - sin 0 
f 2 - 1 - cos 0 
f3 = sin 0 
f 4 =■ cos 0 

fx = sinh 0-0 
f 2 a cosh 0-1 
f 3 = sinh 0 
f 4 = cosh 0 


(A8) 


(A9) 


Here 0 is the change in eccentric anomaly (E - E 0 ) or in hyperbolic anomaly 
(F - F 0 ) , whichever is appropriate. 


A series is used for evaluating f x and f 2 j then f 3 and f 4 can be com- 
puted by use of the defining equations. The series for fx and f 2 are given 


by 


fl = 


n+i 


0 


2n+i 


n = 1 


5r) 


(2n + 1) ! 


(A10) 



m=i 


The factor (-a/|a|) automatically changes the series from circular to hyper- 
bolic functions. This formulation is similar to the one used in reference 10 
(Herrick's variable), but is not valid for the parabolic case. 
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If equation (A7) is to be solved for the value of 0 corresponding to a 
given At, it is necessary to use an iteration. The iteration method used in 
this study can be understood best by examination of the listing of SUBROUTINE 
LAPLCE in appendix D and will not be discussed here. 


The transition matrices were obtained in closed form by differentiating 
equations (Al) with respect to the initial components of position and velocity. 
For convenience, the following notation is defined: 


Sr Sr 
SR q ' SV G 

SV SV 
SRq ' SV D . 


3x3_ matrices of first partial derivatives of £omponents of 
R or V with respect to the components of R Q and V Q 


gradient operator with the components of position as the inde- 
pendent variable 


gradient operator with the components of velocity as the inde- 
pendent variable 


Using these definitions, we can write the transition matrix in partitioned 
form as 




SR 

Sr 

*1 ^ 


SR 0 

Sv 0 

9 3 <P 4 


5v 

Sv 


5R 0 

Sv 0 


From equation (Al) 

cp x = f I + R D (v R f) T + V G (v R g) T 

cp 2 = g I + R Q (V v f) T + VoCVyg) 1 ’ 

9 3 = f I + R D (v R f) T + V G (v R g) T 

cp 4 = g I + R Q (v v f) T + V 0 (v y g) T 

where I is a 3x3 unit matrix. 






(All) 


(A12) 
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it can be seen that 


From the definitions of the fq 

dfi = f 2 d 0 
df 2 — f 3 d 0 

df 3 = f 4 de 

-a 


di’4 = 

Using equations (A13) , one can write 

„2 


f 3 d0 


Vf = 


~ 1 a I 

R^ 


f 3 V0 + 


R o 2 


Vg = - 


Vf = - 


7 


fpV0 - fiVl 


J 


a 


j, 


a 


OT f4V0 + ~R^R f3VR ° + TTT V^ |a 


J, 


M- 1 a | 


s/q 


R, 


vg = g- Vf - 


1 - f 


KE„ <®o ^ (1 - f)(^)K 


ft \ ■n. 


R 1 


(A 13 ) 


(Aik) 

(A15) 


^ (vr) t r (ai6) 


(A1'7) 


The notation VR is used here to mean c®/ OR 0 or BR/5 V q whichever is 
appropriate* Finally, if At is held constant, then 


v(At) = 0 
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This expression can be solved for VQ to give 

AtV 


VO = 


h 


M- 


M- 


R 


fsV 


+ fpV I 


W 


M 


(A18) 


For convenience in programming, two one -dimensional arrays labeled Q and C 
are defined. The part of the Q array which depends only on the values of 
R 0 and V 0 (i.e., not on At or 0) is computed in SUBROUTINE CNSTNT while the 
remainder is computed in SUBROUTINE LA.PLCE. These arrays are defined below in 
tabular form. Note that the elements are given in numerical order and not in 
the order in which they are computed in program. 


Index 


Q array 


0 array 


^15(^3 ~ Q22C17Q2) + Ql7Q3 2 /Q9Ql€ 


2 


Q 22 Q 7 C 17 Q 2 

3 

*3 

Ql7(Q3Ci8 “ 20x70 g) 

k 

U 

^15^6 “ Q 4 Q 7 /Q 9 Q 16 “ Q 7 /O 

5 

f 

Q 2 Q 4 Q 1 1 /Q 9 Q 1 5 

6 

g 

Q 4 C 18 / Q 9 Q 13 Q 16 + Q 7 C 17 

7 

f 

-Q 7 / 0l5 

8 

♦ 

g 

CisCxo + c P 

9 

R o 

Qi 4 Qi 1 Qi 3 Q 2 2 / Qi 6 

10 

a 

Ql4( 02^X8 “ 3QlCi7) 

ll 

|a| 

CiCis - (l - Qs)/QgQi6 

12 

M- 

^2^16 

13 

l/s/l 1 1 a 1 

C 3 Ci6 

ik 

4 1 a l 3 /m- 

R q (1 - f )/R 3 = (1 - 0a)/Q: 

15 

R2 

P-/ R o 3 - Q 12 /O 9 O 21 

16 

R 

R q / r ~ O 9 /Q 16 

17 

la|/R 0 

a /p- = O 10 /O 12 

18 

R 0 /|a| 

Ql7^16^17( 3Ql + Qx8Q3 + 2( 

19 

(Ro •V 0 )/ N /|i a| 


20 

<1 

O 

I\) 


21 

R o 2 


22 

-a/|a| 


23 

At 
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It can be shown that formulas for the gradients in equations (Al4) through 
(Al8) which apply for both elliptic and hyperbolic orbits are as follows: 



v vCr) ■ RT 7 ° - " 2QiaCl7l?0 



®sAl a l p Qio - 

R| ° QisQsQ^i ° 



- QioQiiQisVq 





- 3Qi4C 17 V 0 



a(R 0 • V 0 ) — 

r 3 0 7 ^t 


Qi3V 0 - QigCi5C 1 7R 0 



a(R p - V 0 ) 

M- v/m- I a l 


0l3R o - Qi9CiyV 0 


2k 




V 





a 



f 2 + 


a a 



B 


o 


M* I a | 


[a] 

R 


f 2 V r 


-(« 


Cx5Gi8 + 


03 \ — 


Q 9 Q 1 S, 


*0 " 


Q 11 Q 13 Q 2 — 


Qie 


V r 


VR 


e 




Rq ♦ Vq 

JT Taf 


f 2 v 0 


1 |a| 

VRT R 



v 


O 


— Q 11 Q 3 Q 2 _ 

• ' ClsTc " ~~ aTT" R ° 


Substitution of the gradients into equations (Aik) through (A17) gives 


V-gf = “Qi7Q3Vp0 


^QitCisCiyQsRq + 


Q 17 Q 2 s* 
Qai 0 


= CiR 0 + 


V “ "^IvQ 3VyO + 2Q 17 2 Q 2 Q 1 aCi7V 0 = C 2 R Q + ^3^0 


^Rg = “Qi 4Q,2V^0 - 3QiQi4G 15 C 17 R 0 = CqR g + C 9 V 0 


Vyg = -Ql4Q2Vy9 - 3Qi4Ci 7 QiV q = GqRq + CxoVq 
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V ■ QisfeQiS ' S' B ° + ,7Cis0itRo + ° 7 (i^ B 


_ T_ 

= C 4 R 0 + C 5 v 0 + C t 9i R 


vJ- = Vy0 + Q 7 Ci 7 V 0 + C 7 9 2 T R = CsRq + C S V 0 + Cy^R 

v Q 13 Q 9 Q 1 S 


. ( 1 - Qn) _ T— — — T_ 

V R g = C 16 V R f - - R 0 + C l4 9i R = C n R 0 + C 12 V 0 + Ci^i R 

Q9QlS 


T— — — T— 

Vyg = CisVyf + Ci49 2 R = Ci^Rq + Ci 3 V 0 + Cl 4^2 R 

Finally, substitution in equations (A2) gives 

cp x = fl + R D (V R f) T + V 0 (7 R g) T 

= fi + r 0 (CiR 0 T + c 2 v q T ) + V 0 (CsR 0 T + C 9 V 0 T ) 

= fl + (CxRo + C 8 V 0 )R 0 T + (02 R o + C 9 V 0 )V 0 T 

9 2 = gl + R 0 (V v f) T + V 0 (Vyg) T 

= gl + (c^ 0 + C 9 V 0 )R 0 T + (C 3 R 0 + c 10 v 0 )v 0 T 

9 S = fi + R 0 (v R f) T + Vo(V R g) T 

= 1*1+ C4RQR0 + Cs^o^o + CuVqRo + C 12 V 0 V 0 
+ CyR c ^ T 9 1 + C^V^V 

= fi + (C 4 R 0 + ChVo)R 0 T + (C 5 R 0 + C 12 V 0 )V q T 
+(C 7 R q + C l4 Vo)R 9 i 
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> 

I 

q> 4 = gi + r q ( V) T + v 0 (v v g) T 

= gl + C5R0R0 + CeRoVo + C12V0R0 + C3.3V0V0 
+ CyRjMPs + CuyjFVs 

= gl + (C 5 R 0 + C 12 V 0 )R 0 T + (C S R 0 + C 13 V 0 )V 0 T + (CyR 0 + C l 4 V 0 )R T cp 2 
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APPENDIX B 


EQUATIONS FOR PERTURBING ACCELERATIONS 


The perturbing accelerations are computed on the basis of the four -body 
equations of motion given in reference 7* The coordinate system is Cartesian 
with the X axis positive toward the vernal equinox, the Z axis parallel to 
the earth's polar axis, and the Y axis completing a right handed orthogonal 
set- When the vehicle is within 66,000 km of the moon a selenocentric system 
is used while otherwise it is geocentric. The vehicle's equations of motion 
in the geocentric system are as follows: 


X 



x M ) 


l a m^m U S (X - X g ) p s X s 


Y = 



-x m ) 


Hm Y m 

R 3 
il m 


U S (Y - Y s) 


^s Y s 


Z 





A 3 
^m 


Um^m p s (Z - Z g ) P S X S 
R m 3 " Ag 3 ' R s 3 


where (Xm, Ym, Zm) and (X s , Y s , Z s ) are the positions of the moon and sun, 
respectively- Other quantities are defined as follows: 


R e = Jx 2 + Y 2 + Z 2 
~ s/x m 2 + Y m 2 + Z m 2 
R s = Jx g 2 + Y s 2 + Zg 2 
An =7(X - X m ) 2 + (Y - Y m ) 2 + (Z - Z m ) 2 
Ag =J(X - X s ) 2 + (Y - Y s ) 2 + (Z - Z s ) 2 

p e = 3.98603IXIO 5 km^/sec 2 
Urn = 4 . 8938269XIO 3 km 3 /sec 2 
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F 


p s = 1.3255X10 11 km 3 /seG 2 

<JRq 2 =. 6.6o1j-3958x10 4 km 2 


where Rq is the mean equatorial radius of the earth. 

The three components of the vector u of perturbing accelerations 
(eq. (2)) are obtained by subtracting the first harmonic term (the two-body 
portion) of the acceleration from the total. Thus, in geocentric coordinates: 


U4 =■ X + 


u 5 = Y + 


Ug = Z + 


T? 3 

n e 

p e Y 

Re 3 

p e z 

Re 3 


In selenocentric coordinates , the terms involving J are omitted and 
X, Y, and Z now denote the vehicle's position with respect to the center of 
the moon. In this case 


1x4 


- X + 


A m 3 


u 5 



+ 


^rrr^ 

A m 3 


~ Z + 


^m z 

Am 3 


(Note that in the program listings the position vectors of the moon PM 
and the sun PS and the velocity vector of the moon VM always refer to 
geocentric coordinates . ) 
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APPENDIX C 


SIMPSON'S RULE FOR UNEQUAL TIME INCREMENTS 


It is desired to integrate the transformed forcing function y over the 
time interval from t Q to t 2 . We assume that each component y^ can be fit 
over the interval with an equation of the form 


y = a 2 (t - t Q ) 2 + a x (t - t Q ) + a Q 


(Cl) 


The subscript , i, indicating the component of y has been omitted and equa- 
tion (Cl) may be taken as representing any of the six components of y. 

Let t x be some time between t 0 and t 2 and define 


h =• t - t 0 

hi = ti - t Q 

h 2 = t 2 - t Q 

y 0 = y(t D ) 

yi = y(ti) 


(C2) 


y 2 = y(t 2 ) J 

then equations (Cl) and (C2) can be combined to give 


y 0 “ a o 

yi = a 2 hi 2 + a! hi + a Q \ 
y 2 = a 2 h 2 + a x h 2 + a Q 


(C3) 


30 




+ y 0 h 2 


Collecting terms in y changes equation (c6) to 


' t 2 h 2 hi(2h 2 - 3hi)y 2 + h 2 3 y! + h 2 (3hi - h 2 )(h 2 - hi)y Q ( . 

y dt = ■ - — — ^ l ) 

b Q 6hi(h 2 - hi) 


The coefficients of the y n are functions only of time and are valid for each 
component of y. They are stored in the AI array in the program. 
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APPENDIX D 
PROGRAM LISTINGS 


The following pages present the Fortran IV listings of the main program 
and subroutines used for Danby f s integration method. The single precision 
version is presented here, hut the program is set up for easy conversion to 
double precision computation of the oscillating conic, to more than one inte- 
gration step between rectifications, and for inclusion in a complete guidance 
and navigation program. For this reason all common data are stored in 
"labeled common," and storage is provided for a number of superfluous vari- 
ables. Comment cards have been inserted at appropriate points to explain the 
operation of the program. 

In addition to the subroutines for which complete listings are provided, 
the program uses five general purpose subroutines. Only enough listing is 
presented to explain the use of these subroutines. 
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$ I BFTC TG0503 NOREF 

CTG0503 DANBY INTEGRATION PROGRAM J D MCLEAN 

C OSES SOBROUTINES TG050B, D, S, T, 0, V, W, X, AND Y 

DIMENSION GP(6) ,ED(6) ,DELR{ 3) ,DELV(3 ) ,GC(6,3 ) ,DELX(6) ,XPREV(6) , 
1DELYP(6),DD(6) ,AI (3) ,DELY(6) ,H(3) , AP ( 6 ,6 ) , X ( 6 ) , XPRES(6) , 

2XZER0(6),KK(6),VARLST( 14) ,TX(3) , RR ( 3 ) , VV ( 3 ) , AE < 6 ,6 ) 

COMMON /DNBY/ RZERO( 3 ) , VZERO ( 3 ) ,PM(3),VM(3),Q(23),PS(3),G(3),K(6), 
1ACCREC,ACCUI ,ACCLI , US , UE, UM, C J2 ,DELM ,R , V ,DELS ,DM ( 3 ) , TMSTOP , 
2TIMEA,TIME,CHOVER,HT,HE,RPRES(3) ,VPRES(3) , ESTOP , DELTE , AA ( 3 , 3 , 5) , 

3 1 PRINT, A0(6,6) , I CYCLE, I MOON, I SON , TMOON ( 32 , 2 5 ) , TSUN (4 , 16 ) , DE ( 3 ) 

EQO I VALENCE (XZEROt 1),RZER0( 1) ) , ( XZERO (4 ) , VZERO (1) ) , ( DELX ( 1 ) , DELR 
1(1) ), (DELX(4),DELV( 1) ) , ( XPRES ( 1 ) , RPRES ( 1 ) ) , ( XPRE S (4 ) , VPRES ( 1 ) ) 

C ACCOMULATE TOTAL INTEGRATION TIME IN DOOBLE PRECISION 

DOUBLE PRECISION TAPRS, TADP 

DEFINITION OF VARIABLES 
TS=STOP TIME 

TA=TIME FROM INJECTION TO PRESENT 
UE , UM, US = MU OF EARTH, SUN, MOON 

HT=+1.0 FOR FORWARD INTEGRATION, -1.0 FOR BACKWARD 
HE= STARTING INCREMENT IN ECCENTRIC ANOMALY 
RR= POSITION VECTOR 
VV= VELOCITY VECTOR 

K ( 1 ) = — 1 FOR MOON CENTERED COORDINATES, +1 FOR EARTH CENTERED 
K ( 2 ) NOT USED 

K ( 3 ) = LIMIT ON NUMBER OF PASSES THROUGH INTEGRATION LOOP 
K ( 4 ) = ZERO EXCEPT WHEN COORDINATES ARE TO BE SWITCHED 
K ( 5 ) = ZERO EXCEPT WHEN TRAJECTORY IS COMPLETED 
K ( 6 ) NOT USED 

DATE= JULIAN DATE OF START OF SUN-MOON TABLES 
DAYS= NOT USED 
ACCREC= A SUB R -- SEE TEXT 
ACCUI = A SUB U 
ACCL I = A SUB L 

C J2= J*(R SUB Q)**2 — SEE APPENDIX B 
CHOVER = RANGE FROM MOON OF COORDINATE SWITCH 
T I = INJECTION TIME 
PS= POSITION VECTOR OF SUN 
PM= POSITION VECTOR OF MOON 
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DELTE IS TOTAL INCREMENT IN ECCENTRIC ANOMALY FROM LAST 
RECTIFICATION 

VM=VELOC I TY VECTOR OF MOON 
DELM= VEHICLE-MOON DISTANCE 

AE IS TRANSITION MATRIX FROM START OF TRAJECTORY TO PRESENT 
CALL CLOCK { TX ) 

C FROM HERE TO STATEMENT 100 INITIALIZES PROGRAM FOR ONE TRAJECTORY 

C INTEGRATION 

200 READ! 5, 111) TS , TA , UE I , UMI , HT I , HE I , ( RR ( I ) , I = 1 , 3 ) , ( VV { I ) , I = 1 , 3 ) 

111 FORMAT ( 3E16.8 ) 

RE AD ( 5 ,90 1 ) (K(J)»J=1,6) 

READ( 5, 111 ) DATE, DAYS » T IME, ACCREC »ACCUI»ACCLI,CJ2,US,C HOVER 

901 FORMAT (615) 

TMSTOP=TS 

TIMEI=TIME-TA 

TIMEA=TA 

T APRS=TA 
TADP=TA 
UE=UE I 
UM=UMI 
HT=HT I 
HE=HE I 

DO 810 1=1,3 
RZERO ( I ) =RR { I ) 

810 VZERO ( I ) = VV ( I ) 

WRITE(6,902)TMSTOP,TIMEA,UE 

902 FORMAT ( 1H 1 , 7HTMST0P= , D25. 16, 6HT IMEA= , D2 5 . 16 , 3HHE= ,D25 . 16 ) 

WRI TE ( 6 , 903 ) UM, HT, HE 

903 FORMAT ( 1H , 3HUM=, 025*16, 3HHT=,D2 5.16, 3HHE=,D25. 16) 

WRITE ( 6 , 904 ) ( RZERO I I ), 1 = 1,3) 

904 FORMAT ( 1H , 2HX=, D25. 16, 2HY= , D2 5 . 16 , 2HZ= , D25 . 16 ) 

WRITE (6,905) (VZERO ( I ), 1=1,3) 

905 FORMAT ( 1H , 3HXD=» D2 5.16,3 HYD=,D2 5.16»3HZD=» 025*16) 

WRITE(6,906) ( K( I ) ,1=1,6) 



906 FORMAT* 1HO , 2HK= , 6 1 5 ) 

WRITE(6, 118) 

118 FORMAT* 1H-,5X,5HTIMEA,12X,5HX(KM) ,12X,5HY(KM) ,12X,5HZ(KM) , 10X,10HX 
ID* KM/SEC ) ,7X, 10HYD(KM/SEC) ,7X, 10HZD* KM/SEC ) ,5X,9HSTEP SIZE) 

WRITE (6, 119 ) 

119 FORMAT* 1H ,6X,4HDELX, 13X, 4HDEL Y , 1 3X , 4HDEL Z , 13X ,4HDELR , 14X , 1HR , 1 6X , 
11HV, 15X,4HDELM) 

CALL MSLOAD(TMDON( 1, 1) , 32,TSUN( 1,1) ,4,DATE,N0MPTS,N0SPTS ) 

CALL CLOCK ( TX ( 2 ) ) 

I M00N=0 
I SUN=0 
KERR=0 
HE P=HE 
I CYCLE=0 
DO 10 1=1,36 

10 AE(I,1)=0.0 
DO 11 1=1,6 
DELX* I )=0.0 

11 AE ( I , I )=1.0 

C BEGINNING OF INTEGRATION LOOP 

100 DO 7 1=1,3 

7 RPRES* I ) =RZERO ( I ) 

CALL MOON* TIME, PM* 1) ,TMOON{ 1, 1 ) ,32,N0MPTS»KERR»3HP0S, I MOON ) 

9 CALL SUNP (TIME, PS* 1),TSUN* 1,1) ,4»N0SPTS,KERR»ISUN) 

8 CALL GVEC 

C GVEC COMPUTES FORCING FUNCTION 

DELMP=DELM 

V=0.0 

DO 6 1=1,3 
6 V=V+VZERO ( I )**2 
V=SQRT ( V ) 

C OUTPUTS TIME IN DAYS 

TIMEB=TI MEA /86400 • 0 

10 3 FORMAT* 1H0 , E 16.8, 6E 17. 8 , E13.4/7E 17.8 , 1 5 ) 

IF ( K ( 5 ) ) 221,221,200 


<J0 

VJ1 


UG 

ON 


C CNSTNT COMPUTES PART OF Q ARRAY AFTER RECTIFICATION 

221 CALL CNSTNT 


DO 21 1=1,3 
21 GC ( I , 1 ) =0 .0 

DELTE=0.0 
ESTOP=l *0 
TPRS=TI ME 
TAPRS=TADP 
212 MA= 1 

12 ICYCLE=ICYCLE+1 

IF ( ICYCLE-KI 3) ) 110,200,200 
110 GO TO (24, 17,17) ,MA 

24 K ( 5 ) =0 

C EACH TRANSFORMED FORCING FUNCTION, Y IN TEXT, IS ONE COLUMN OF GC 

C STORE INITIAL VALUE OF G IN FIRST COLUMN OF GC 

DO 25 1=1,3 

25 GC ( 1+3, 1 )=G { I ) 

241 MA=2 

GO TO 12 

C INCREMENT DELTE ONE STEP 

17 DELTE=DELTE+HE 

C COMPUTE NEW 2 BODY POSITION, RPRES , AND VELOCITY, VPRES 
CALL LAPLCE 


C COMPUTE TRANSITION MATRIX, AO, FROM RECTIFICATION TO PRESENT 

CALL PARTLS 

H { MA-l ) =Q ( 23 ) 

TAPRS=TADP+Q(23) 

TPRS=TAPRS+TIMEI 

CALL SUNP(TPRS,PS( 1) ,TSUN( 1,1) ,4,N0SPTS,KERR,ISUN) 

19 CALL MOON(TPRS,PM( 1) ,TMOON( 1, 1) »32»N0MPTS»KERR,3HP0S, I MOON ) 

C COMPUTE NEW FORCING FUNCTION 

20 CALL GVEC 



C TRANSFORM G TO LAST RECTIFICATION TIME AND STORE IN GC 

CALL XPYXM(G(l),AA(l,l,n,GC(4,MA) ,1,3, 3, 3 , 1 ) 

DO 22 1=1,3 
22 G ( I )=-G ( I ) 

CALL XPYXM( G ( 1)»AA( 1, 1,2) ,GC( 1,MA) ,1,3, 3,3,1) 

C TEST WHETHER GC ARRAY HAS BEEN FILLED 
GO T0(200,230,13) ,MA 

230 MA=3 

GO TO 12 

13 IF ( K { 5 ) »GT *0 ) GO TO 30 

C TEST WHETHER STOP TIME HAS BEEN EXCEEDED, IF SO 
213 TST=HT* ( TMSTOP-TAPRS ) 

IF(TST) 14,18,30 

C COMPUTE EXACT VALUE OF HE TO STOP AT CORRECT TIME 

14 Q(23)=TMST0P-TIMEA 
EST0P=0 .0 
DELTE=0.0 

CALL LAPLCE 
HE=(EST0P-DELTE)/2.0 
K ( 5 ) = 1 
GO TO 241 
18 K ( 5 ) = 1 

C COMPUTE AI ARRAY — SEE APPENDIX C 
30 SS=H(2)-H(1) 

AI ( 1 ) = SS*H(2)*(3.0*H(1)-H(2) ) 

AI ( 2 )=H ( 2)*#3 

AI (3)*H< 1)*H(2)*(2#0*H(2)-3.0*H(1) ) 

SS=SS*H( 1)*6.0 
DO 29 1=1,3 
29 AI ( I )=AI ( I ) /SS 
SS=H(2)/2.0 





C GP= INTEGRAL OF TRANSFORMED FORCING FUNCTION 

CALL XPYXM(GC( 1, 1),AI ( 1 ) , GP ( 1 ) ,6 , 3 ,3 , 1 , 1 ) 

C COMPUTE DG — DG= A SUB Q IN TEXT 

DO 229 1=1,6 

ED (I )=GP(I )-(GC( I,1)+GC( 1,3) )*SS 
229 DD ( I ) =GP ( I )-H(2)#GC(I,l) 

CALL ROOT ( DD( 1 ) , DR2, DDM, 1 } 

CALL ROOT ( DD ( 1),DV2,DDM,4) 

DG=0.0 
DO 40 1=1,3 

40 DG=DG+( ED{ I )**2 ) /DR2+( ED( 1+3 ) **2 ) /DV2 
TST5=DG-ACCUI 

TST6=DG-ACCLI 
IF(TST5.GT.0.0) GO TO 41 
43 DO 243 1=1,6 

C DELX IS USED ONLY WHEN NOT RECTIFYING EVERY STEP 

243 GP( I ) =GP ( I )+DELX( I ) 

C DELY= PERTURBATION STATE VECTOR AT PRESENT 

CALL XPYXM(A0»GP,DELY,6,6,6,1,1) 

C CHECK CONSTRAINTS ON MAGNITUDE OF PERTURBATION 
CALL ROOT (DELY(1),DRSQ,RDEL,1) 

449 IF (R-DELM) 450,450,451 

450 TST7=RDEL/R-ACCREC 
TST8=RDEL/R-0.2*ACCREC 
GO TO 452 

451 TST7=RDEL/DELM-ACCREC 
TST8=RDEL/DELM-0.2*ACCREC 

452 IF(TST7.LE.0.0) GO TO 247 

C THIS SECTION USED WHEN STEP SIZE MUST BE REDUCED 

41 HE=HE/2.0 
251 K ( 5 ) =0 

DELTE=0.0 
GO TO 241 



c 


UPDATE TIME AND AE AFTER RECTIFICATION 
250 K (6 ) = 1 

DO 51 1=1,6 
51 DELX ( I ) =0 .0 
TIME=TPRS 
TADP=TAPRS 
TIMEA=TAPRS 

60 CALL XPYXM(A0,AE,AA,6,6,6,6, 1) 

DO 61 1=1,36 

61 AE ( I , 1 > = AA C 1,1,1) 

C RETURN TO BEGINNING OF LOOP UNLESS TRAJECTORY IS TO BE TERMINATED 

62 IF ( K( 5 ) • LE.O ) GO TO 100 

CALL CLOCK ( TX ( 3 ) ) 

WRITE(6,2) <TX( I ) ,1=1,3) , ICYCLE 

2 FORMAT! 1HQ, 17HEND OF TRAJECTORY ,5X ,2HTX, 3E 17.8 , 5X,6HICYCLE , I 5 ) 

WRITE (6»103)TIMEB»(XZER0(J)»J=1»6)»HE»(DELY(J),J=1»3)»RDEL,R»V,DE 
1LM, ICYCLE 
GO TO 100 


C CHECK FOR COORDINATE SWITCH AND ADJUST HE SO THAT DELM WILL LIE 

C WITHIN 5000 KM OF CHOVER 

247 TST2=K( 1 ) 

TST2=TST2*( DELM- CHOVER ) 

K ( 6 ) =0 

TST3=ABS(TST2)~ 5000.0 
IF ( TST2. LT .0.0) GO TO 39 
IF ( TST3.GT.0.0 ) GO TO 37 
33 SHFT=K( 1 ) 

GO TO 38 

37 RATIO=ABS(DELMP-CHOVER)/ABS( DELMP-DELM ) 

HE=HE*RATIO 

GO TO 251 

38 K(4)=l 
GO TO 47 


LO 

VD 



§ C CHECK WHETHER STEP SIZE SHOULD INCREASE 

39 IF(TST6.GT.0.0) GO TO 47 
IF(TST8.GT.0.0) GO TO 47 
45 HE=HE*2.0 


C THIS SECTION RECTIFIES CONIC AND SWITCHES COORDINATES IF NEEDED 

47 DO 48 1=1,6 

48 XZEROU )=XPRES( I )+DELY( I ) 

DELMP=DELM 

IF ( K ( 4 ) ) 250,250,52 

52 CALL MOON( TPRS, VM( 1 > , TMOONt 1,1) »32»N0MPTS»KERR»3HVEL, I MOON ) 
IF(KERR.NE.O) GO TO 200 

53 DO 54 1=1,3 

XZERO ( I )=XZERO(I )+SHFT#PM{ I ) 

54 XZERO ( I+3)=XZER0( I+3)+SHFT*VM( I ) 

K ( 4 ) =0 

IF(SHFT.GT.O.O) GO TO 56 

55 K ( 1 ) = 1 
WRITE (6,102) 

102 FORMAT! 1H0,59HTHE FOLLOWING SETS OF DATA ARE IN MOON CENTERED COOR 
1DINATES ) 

GO TO 57 

56 K( 1 ) = — 1 
WRITE (6,101) 

101 FORMAT ( 1HO,60HTHE FOLLOWING SETS OF DATA ARE IN EARTH CENTERED COO 
1RDINATES ) 

57 HE=HEP 

GO TO 250 
END 
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$ I BFTC TG050V NOREF 

SUBROUTINE LAPLCE J D MCLEAN COMPUTES TWO BODY POSITION AND 
VELOCITY 

SUBROUTINE LAPLCE 

COMMON /DNBY/ RZ ERO{ 3 ) , VZERO ( 3 ) , PM ( 3 ) , VM { 3 ) , Q ( 23 ) , PS ( 3 ) , G ( 3 ) , K ( 6 ) , 
1ACCREC,ACCUI ,ACCLI,US,UE,UM,CJ2,DELM,R,V,DELS,DM(3) ,TMSTOP, 
2TIMEA,TIME,CHOVER,HT,HE,RPRES(3) ,VPRES (3) ,ESTOP, DELTE,AA{ 3,3,5) , 
3IPRINT,A0(6,6),ICYCLE,IM00N,ISUN,TM00N(32,25),TSUN(4,16),DE(3) 

E IS INCREMENT IN ECCENTRIC ANOMALY 
E=DELTE 
EP=DELTE 
DTP=0.0 
EINC=HE 


DO 15 1=1,20 

IF ESTOP IS ZERO KEPLERS EQUATION IS SOLVED TO GET DELTE FOR STOP 
AT CORRECT TIME, OTHERWISE OUTER LOOP IS OMITTED 
I F ( ESTOP ) 9,8,9 

8 E=EP+E INC 

C KD=0 UNLESS EINC IS TO BE SUBTRACTED FROM E 
KD=0 

10 IF(E-EP) 9,14,9 

9 Q(1)=0.0 
Q ( 2 ) =0*0 

C S IS A DUMMY VARIABLE 
S=1.0 


•p- 

H 


C THIS LOOP SOLVES SERIES FOR Q( 1 ) AND Q(2) — SEE APPENDIX A 
DO 11 J=l, 50 
RN= J 

RN=2«0*RN 

S= S*Q(22)*E**2/(RN*(RN-1.0) ) 

H=Q(22)*S 
QZ=Q ( 2 ) 

Q(2)=Q(2)+H 

41 Q ( 1 )=Q ( 1 ) +H*E/( RN+1.0 ) 



on on 


-£■ 

ro 


C 


SERIES IS CARRIED OUT UNTIL LAST TERM DOES NOT CHANGE Q ( 2 ) 

I F ( Q ( 2 ) — Q Z ) 11,12,11 
11 CONTINUE 

WRITE (6, 101)S 

101 FORMAT ( 2H0 , 38HQ 1-Q2 SERIES HAS FAILED TO CONVERGE S=,E15.8) 


C COMPUTE TIME INCREMENT 

12 Q(4)=1.0+Q(2)*Q(22) 

Q(3)=E+Q( 1 ) *Q ( 22 ) 

42 S=Q< 14)*<Q( 1)+Q( 18)*Q(3)+Q( 19)*Q(2) ) 

C IF NOT ITERATING FOR STOP TIME LEAVE LOOP 

I F { ESTOP ) 14,23,14 

C TEST WHETHER TIME INCREMENT HAS BEEN FOUND TO 7 PLACES, IF NOT 

23 TST 1= ABS( 1.0-S/QI23) J-.1E-6 
IF(TSTl) 14,14,13 

C TEST WHETHER TIME INCREMENT IS TOO LARGE, IF NOT 

13 TST2= ABS ( S ) -ABS ( Q ( 23 ) ) 

IFITST2) 30 ,30,31 

C INCREMENT E AND REPEAT ITERATION 

30 EP=E 
DTP= ABS ( S ) 

GO TO 15 

IF E IS TOO LARGE AND KD=0 COMPUTE RATIO OF ACTUAL TIME INCREMENT 
TO DESIRED ONE 

31 IF(KD) 34,34,38 

34 RATI 0= ( ABS ( Q ( 23 ) ) -DTP ) / ( ABS(S)-DTP) 

IF RATIO IS LESS THAN 0.1 MULTIPLY EINC BY 0.1 AND START AT 
BEGINNING OF LAST INTERVAL 
IF (RATIO-O. 1 ) 35,35,36 
E INC=E INC*0. 1 
GO TO 8 


35 



IF RATIO IS GREATER THAN 0.9 MULTIPLY EINC BY 0.1, SET KD=1 AND 
BEGIN REDUCING E 

36 IF (RAT 1 0-0. 9 ) 39,37,37 

37 E INC=E INC*0. 1 

38 E=E-E INC 
KD= 1 

GO TO 10 


IF RATIO LIES BETWEEN 0.1 AND 0.9 USE LINEAR INTERPOLATION 
39 TINC=TINC*RATIO 
GO TO 8 

15 E=E 

WRITE (6,102)TST1 

102 FORMAT { 1H0»43HTIME ITERATION HAS FAILED TO CONVERGE TST1=,E15.8) 

IF ITERATING FOR STOP TIME RETURN TO MAIN PROGRAM 
14 IF (ESTOP ) 25,24,25 

24 ESTOP=E 
GO TO 33 

COMPUTE NEW POSITION AND VELOCITY 

25 Q ( 23 ) = S 

ROA=Q (2)+Q( 18 ) #Q ( 4) +Q ( 19 ) *Q ( 3) 

Q ( 5 ) = 1.0-Q ( 17)*Q(2) 

Q ( 6 ) =Q ( 23 ) -Q ( 14 ) *Q ( 1 ) 

Q(7)=-Q( 17)*Q(3)/(R0A*Q( 14) ) 

Q ( 8 ) = 1.0-Q ( 2 ) /ROA 
Q ( 16 ) =R0A*Q (11) 

43 Q( 15 )=Q( 16)**2 
100 DO 16 1=1,3 

RPRESU )=Q( 5)*RZER0( I ) +Q ( 6 )*VZERO ( I ) 

16 VPRESU )=Q(7)*RZER0( I )+Q( 8)*VZER0( I ) 

33 RETURN 

END 
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P $1 BFTC TG050W NOREF 

SUBROUTINE CNSTNT J D MCLEAN COMPUTES CONSTANTS FOR LAPLCE 

SUBROUTINE CNSTNT 

COMMON /DNBY/ RZERO ( 3 ) , VZERO ( 3 ) , PM { 3 ) , VM ( 3 ) ,Q ( 23 ) , PS { 3 ) , G ( 3 ) , K ( 6 ) , 
1ACCREC, ACCUI , ACCl I ,US , U.E , UM, CJ2 , DELM ,R , V ,DELS , DM ( 3 ) , TMSTOP , 
2TIMEA,TIME,CHOVER,HT,HE,RPRES(3) ,VPRES (3) , ESTOP , DELTE , AA ( 3, 3, 5 ) , 
3IPRINT,AO(6,6) , ICYCLE, IMOON, ISUN,TMOON(32 ,25 ) , TSUN (4 , 16 ) , DE ( 3 ) 

I F ( K ( 1 ) ) 17,17,18 

17 Q ( 12 ) =UE 
GO TO 19 

18 0(12) =UM 

19 Q(19)=0.0 
Q ( 20 ) =0 *0 
Q ( 21 ) =0.0 
DO 10 1=1,3 

Q( 19)=Q( 19)+RZER0( I )*VZERO( I ) 

Q(20)=Q(20)+VZER0( I )**2 
10 Q ( 21 )=Q ( 21 ) +RZERO( I )**2 
Q ( 9 ) = SORT ( Q ( 21 ) ) 

Q( 1)=C2.0/Q(9) ) — (QC20)/Q( 12) ) 

Q( 10) = 1.0/Q( 1) 

Q( 11)= ABS(Q( 10) ) 

Q ( 1 3 ) = 1 .0 / { SORT (Q( 11)*Q( 12) ) ) 

Q(22)=-Q(10)/Q( 11) 

Q ( 17 )=Q ( 11 ) /Q (9 ) 

Q ( 18 ) = 1.0/Q ( 17) 

Q( 14)=Q( 13)*Q( 11)**2 
40 Q(19)=Q(19)*Q(13) 

RETURN 
END 
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$1 BFTC TG050X NOREF 

C SUBROUTINE GVEC COMPUTES FORCING FUNCTION — SEE APPENDIX B 

SUBROUTINE GVEC 
DIMENSION DS ( 3 ) , D2 ( 3 ) 

COMMON /DNBY/ RZEROt 3) , VZERO( 3) , PM ( 3 ) , VM ( 3 ) ,Q { 23 ) ,PS ( 3 ) ,G { 3 ) , K ( 6) , 
1ACCREC,ACCUI,ACCLI , US, UE , UM, CJ2 ,DELM ,R , V ,DELS ,DM ( 3 ) ,TMSTOP, 
2TIMEA,TIME, CHOVER,HT, HE, RPRES { 3) , VPRES ( 3 ) , ESTOP ,DELTE , AA ( 3 , 3, 5 ) , 

3 1 PRINT, AO (6, 6) , I CYCLE, I MOON, I SUN ,TMOON ( 32 ,25 ) , TSUN (4, 16 ) , DE ( 3 ) 

DM= VEHICLE-MOON POSITION VECTOR 
DE=VEHICLE-EARTH POSITION VECTOR 

CALL ROOT (PM(1),RMSQ,RM,1) 

RM3=RMSQ*RM 
I F ( K ( 1 ) ) 10,10,12 

10 DO 11 1=1,3 
DM( I )=RPRES ( I )-PM( I ) 

11 DE ( I )=RPRES ( I ) 

GO TO 14 

12 DO 13 1=1,3 
DM ( I ) =RPRES ( I ) 

13 DE C I )=RPRES(I ) +PM ( I ) 

14 DO 15 1=1,3 
DS(I)=DE(I )— PS ( I ) 

15 D2( I )=PMt I ) — P S ( I ) 

GALL ROOT (DE(1)»RSQ»R,1) 

CALL ROOT { DM C 1),DELMSQ,DELM,1) 

CALL ROOT ( DS(1), DELS Q, DELS, 1) 

CALL ROOT(D2{l),DEL2SQ»DEL2»l) 

DELS3=DELSQ*DELS 
U3=— US/DELS3 
R3=RSQ*R 

I F { K ( 1 ) ) 16,16,17 

16 DELM3=DELMSQ*DELM 
CALL ROOT C P S ( 1) ,RSSQ,RS»1) 

RS3=RSSQ*RS 
Ul=— UM/DELM3 
U2=-UM/RM3 
U4=-US/RS3 



GO TO 18 

17 DEL23=DEL2SQ*DEL2 
U1=-UE/R3 
U2=UE/RM3 
U4=US/DEL23 

18 DO 19 1=1,3 

19 G(I)=U2*PM(I)+U3*DS(I) 

IF ( K ( 1 ) ) 20 , 20 , 22 

20 BRK=1*0-5.0*RPRES(3)**2/RSQ 
C0=-UE*CJ2/(R3*RSQ) 

DO 21 1=1,3 

21 G(I)=G(I)+DE(I)*BRK*CO +U1*DM ( I ) +U4*PS { I ) 
G(3)=G(3)+2.0*DE(3)*C0 

GO TO 24 

22 DO 23 1=1,3 

23 G ( I )=G(I)+U1*DE(I )+U4*D2(I ) 

24 R=R 
RETURN 
END 



$1 BFTC TG050Y NUREF 

C SUBROUTINE PARTLS J D MCLEAN COMPUTES TWO BODY TRANSITION 

C MATRICES — SEE APPENDIX A 

SUBROUTINE PARTLS 

DIMENSION QR(5) ,C( 18) ,QQ{ 5) ,S(3) , D2 ( 4 ) ,D3(4),D4(4) 

EQUIVALENCE ( D2 { 1 ) , C ( 2 ) ) , ( D3 ( 1 ) , C ( 8 ) ) , ( D4 (1) , C ( 9 )) 

COMMON / DNBY / RZ ERO ( 3 ) , VZE RO ( 3 ) , PM ( 3 ) , VM ( 3 ) , Q ( 23 ) , PS ( 3 ) , G ( 3 ) , M ( 6 ) , 
1ACCREC,ACCUI , ACC LI , US , UE , UM, C J2 , DELM ,R , V ,DELS , DM ( 3 ) ,TMSTOP, 
2TIMEA,TIME,CHOVER,HT,HE,RPRES( 3) , VPRES < 3 ) , ESTOP , DELTE , AA ( 3 , 3 , 5 ) , 

31 PRINT, AO (6, 6) , I CYCLE, I MOON, I SUN , TMOON ( 32 , 25 ) , TSUN ( 4 , 16 ) , DE ( 3 ) 

DO 90 1=1,5 
QR ( I ) =Q ( I +3 ) 

90 QQ ( I ) =Q (I +4) 

C( 17 ) =Q ( 10) /Q ( 12) 

C( 16)=Q(9)/Q( 16) 

C( 15 )=Q( 12) /(Q(9)*Q( 21 ) ) 

C( 14)=1.0-Q(8) 

C(7)=-Q(7)/Q( 15) 

C ( 9 ) =Q ( 17 ) #C ( 16) 

C(2)=Q(2)*C( 17) 

C(5)=Q(4)/(Q(9)*Q( 13)*Q( 16) ) 

C { 10 ) = 3.0*C(17)*Q( 1) 

C( 12)= Q I 3 ) *C ( 17) 

C( 18)*C(9)*(C( 10 )+Q( 18)*C( 12)+2.0*g(19)*C(2) ) 

C ( 3 ) = Q( 17)*(Q(3)*C( 18)-2.0*C( 2) ) 

C( 1)=(C{3)-Q(22)4C(2) )*C(15)+C(9)*Q(3)**2/Q(21) 

C(6)=Q(7)*C( 17)+C( 5)*C( 18) 

C( 10) = -Q( 14)*(C( 10)-U(2)*C( 18) ) 

C(4)=C(6)«C( 15)-C(5)*Q( 13 )*Q( 7)-Q( 7 ) /Q ( 21 ) 

C(9)=-C{ 14)*C(2)*Q(22) 

C ( 2 ) = C ( 2 ) *Q ( 22 ) *Q ( 7 ) 

C(5)=C( 5)*C( 14)*Q( 13) 

DO 9 1=1,3 

9 C ( I + 10 ) =C ( I )*C( 16) 

C(11)=C( 1 1 ) -C ( 14) /Q ( 21 ) 

C(14)=C(14)/Q(15) 

C(8) = C( 15)*C( 10 ) +C ( 2 ) 



-t=- 

OD 


DO 12 K=l, 5 
IFIK-3) 10,12,10 

10 DO 1 1 1 = 1,3 
DO 51 J= 1 , 3 

51 AA( I, J,K)=(C(K)*RZERO( I ) +D3 ( K ) *VZER0 ( I ) )*RZERO( J) + 
1 ( D2 ( K ) *RZ ERO ( I ) + D4 ( K ) * VZERO ( I ) ) *VZERO ( J ) 

IF ( K — 3 J 53,12,52 

52 AA(I»I»K)=AA(I,I»K)+UR(K) 

GO TO 11 

53 AA(I,I,K)=AA(I,I,K)+QQ(K) 

11 CONTINUE 

12 CONTINUE 

DO 41 1=1,3 

41 S ( I )=C( 7)*RZER0( I ) +C( 14)*VZER0( I ) 

DO 13 1=1,3 
SR=0 .0 
SV=0.0 
DO 14 J= 1 , 3 

SR = SR+RPRES( J ) =4 AA ( J, I ,1) 

14 SV=SV+RPRES(J)*AA(J,I,2) 

DO 13 K= 1 , 3 

AA(K,I,4)=AA(K,I ,4)+S(K)4SR 

13 AA(K,I,5)=AA(K,I,5)+S(K)*SV 
DO 17 1=1,3 

DO 17 J=l, 3 

AO ( I , J ) = AA ( I , J, 1 ) 

AO ( 1 + 3, J + 3)=AA( I , J, 5) 

AO ( I , J + 3 ) = AA ( I , J , 2 ) 

17 A0( 1 + 3, J ) = AA ( I , J,4) 

50 RETURN 
END 



c 


GENERAL PURPOSE SUBROUTINES 


$1 BFTC TG050B NODECK, NOREF 
C MATRIX MULTIPLY ROUTINE 

CTG050B SUBROUTINE XPYXM J D MCLEAN J=1 C=A*B J=2 C=A*6T J=3 C=AT*B 
SUBROUTINE XPYXM ( A, B , C, NR A, NCA ,NRB ,NCB , J ) 

C NRA=NO. OF ROWS IN A 

C NCA= NO. OF COLUMNS IN A, ETC 


$ I BFTC TG050D 

CTG050D SUBROUTINE ROOT J D MCLEAN 
SUBROUTINE ROOT ( X, RSQ, R , NFCA ) 

C SUBROUTINE ROOT COMPUTES MAGNITUDE, R, AND SQUARE, RSQ, OF VECTOR 

C GIVEN BY ANY 3 CONSECUTIVE TERMS OF ARRAY, X, STARTING WITH 

C X (NFCA ) 


$1 BFTC TG050S 

C LOADS COEFFICIENTS FOR MOON AND SUN POLYNOMIALS L MCGEE 

SUBROUTINE MSLOADf TMOON, LTABM, TSUN, LTABS , DATE ,KM, KS ) 

C THE CALLING PROGRAM WILL ORDINARILY HAVE TMOON DIMENSIONED AS 

C TMOON( LTABM, 25) AND TSUN DIMENSIONED AS TSUN ( LTABS, 16 ) WHERE 

C LTABM AND LTABS ARE THE NUMBER OF ROWS IN THEIR RESPECTIVE TABLES. 


$ I BFTC TG050T 

C COMPUTES SUN POSITION FROM POLYNOMIAL COEFFICIENTS L MCGEE 

SUBROUTINE SUNP( TIME, PS, TSUN, LTABS , NOS PTS, KERR, I ) 

C PS IS POSITION VECTOR OF SUN, TIME IS FROM START OF TABLE IN SEC. 


$ I BFTC TG050U 

C COMPUTES MOON POSITION AND VELOCITY FROM POLYNOMIAL L MCGEE 

SUBROUTINE MOON (TIME, QM, TMOON, LTABM, NOMPTS, KERR, POV,I ) 

C QM IS POSITION OR VELOCITY VECTOR OF MOON DEPENDING ON WHETHER 

C POV IS POS OR VEL IN CALLING PROGRAM 
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Figure 1.- Total number of integration steps with step size adjusted o nl y by 

constraint on r/R. 
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Figure 2.- Total number of integration steps with step size adjusted by 

constraints on r/R and Aq. 
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Figure 3.- Deviations in X component of position. 




Figure 4.- Accuracy and integration time for transearth trajectory. 
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Figure 5.- Accuracy and integration time for translunar trajectory. 
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